rm(list=ls())
set.seed(08536)
if(!is.null(sessionInfo()$otherPkgs)){
  invisible(lapply(paste0('package:', names(sessionInfo()$otherPkgs)), detach, character.only=TRUE, unload=TRUE))
}
if (!require(pacman)){
  install.packages(pacman)
  library(pacman)
}
p_load(stargazer,
       survey,
       dplyr,
       questionr,
       srvyr,
       jtools,
       interactions,
       ggplot2,
       coefplot,
       corrplot,
       sjPlot,
       interplot,
       ggeffects,
       viridis,
       mice,
       labelled,
       mitools,
       miceadds,
       dotwhisker,
       texreg,
       gtsummary,
       randomForest,
       caret,
       poliscidata)

##########Helper functions for dealing with multiple imputation################
extract.mi<-function(model){
  sink("file")
  s<-summary(model) 
  names<-names(model$coefficients) 
  co<-s$results
  se<-s$se
  pval<-2*pnorm(abs((co/se)),lower.tail = F)
  nimp <- model$nimp
  # rs<-s$r.squared 
  # adj<-s$adj.r.squared 
  # n<-nobs(model) 
  gof<-c(nimp) 
  gof.names<-c("Num.\\imputations") 
  tr<-createTexreg(
    coef.names=names, 
    coef=co, 
    se=se, 
    pvalues=pval, 
    gof.names=gof.names, 
    gof=gof) 
  sink()
  return(tr)
}
setMethod("extract",signature=className("MIresult","mitools"), definition=extract.mi)

##Read data in
dat <- readRDS("../Data/clean_data.RDS") %>% remove_labels()
dat$hist[is.na(dat$hist)] <- 99
#Not run
# factors <- dat %>% dplyr::select_if(is.factor)
# dat$territ_loss_concern_scaled <- as.vector(dat$territ_loss_concern_scaled)
# datimpute <- dat %>% dplyr::select_if(is.numeric) %>% mice(.,method="cart")
# datimpute <- as.mids(complete(datimpute,"long",include=T) %>% bind_cols(bind_rows(replicate(6,factors,simplify = F)))
# saveRDS(datimpute,"data_imputed.RDS")
svyimpute <- imputationList(complete(readRDS("../Data/data_imputed.RDS"),"all"))
dat4 <- readRDS("../Data/clean_data_noV4.RDS") %>% mutate(territ_loss_concern_scaled = scales::rescale(territ_loss_concern))

p13 <-readRDS("../Data/panel13.rds")
p34 <-readRDS("../Data/panel34.rds")

meplot <- function(model,var1,var2,ci=.95
                   ){
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov(model)
  z0 <- seq(0,100,length.out=101)
  dy.dx <- beta.hat[var1] + 
    beta.hat[(length(beta.hat)-1)]*z0 + 
    beta.hat[(length(beta.hat))]*z0^2
  se.dy.dx <- sqrt(cov[var1,var1] + 
                     z0^2*cov[(nrow(cov)-1),(ncol(cov)-1)] + 
                     z0^4*cov[nrow(cov),ncol(cov)] + 
                     2*z0*cov[var1,(ncol(cov)-1)] +
                     2*(z0^2)*cov[var1,ncol(cov)] +
                     2*(z0^3)*cov[(nrow(cov)-1),ncol(cov)]
                   )
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  return(data.frame(x=z0,est=dy.dx,ci.upper=upr,ci.lower=lwr))
}
makeform <- function(country,outcome,controls=T,interact = F, obs=F,diaspora=F,fe=T,scale=F,quadratic=F){
  #Edit to change baseline controls
  treat <- case_when(diaspora ~ "D",
                     obs&(!scale) ~ "territ_loss_concern",
                     obs&scale ~ "territ_loss_concern_scaled",
                     TRUE ~ "territorial")
  if(interact){
    if(!quadratic){
      basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"), 
                   "age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc", 
                   "male", 
                   "city",
                   "town",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "natid_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01")
    } else{
      basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"),
                   paste(ifelse(scale,"I(territ_loss_concern_scaled^2)","I(territ_loss_concern^2)"),treat,sep="*"),
                   "age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc", 
                   "male", 
                   "city",
                   "town",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "natid_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01")
    }
    
  } else{
    basespec <-c(treat,
                 "age", 
                 "educ_univ_complete", 
                 "educ_sec_univprep", 
                 "educ_sec_techvoc", 
                 "male", 
                 "city",
                 "town",
                 "covid_civlib_restrict_ind",
                 "econ_intervention_ind",
                 "pro_global_ind",
                 "democ_ind",
                 "natid_ind",
                 "neigh_lgbt",
                 "relig_import01",
                 "pol_int01")
  }
  if(!controls){
    if(interact){
      if(!quadratic){
        basespec <- paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*")
      } else{
        basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"),
                     paste(ifelse(scale,"I(territ_loss_concern_scaled^2)","I(territ_loss_concern^2)"),treat,sep="*")
                     )
      }
    } else{
      basespec <- treat
    }
  }
  #Add country-appropriate region fixed effects
  
    geo <- case_when(country == "Romania" ~ "County_Rom",
                     country == "Hungary" ~ "County_Hung",
                     country == "Turkey" ~ "Province_Turk",
                     country == "Germany" ~ "Bundesland_Germ",
                     country == "All" ~ "Country"
    )
  
  
  
  if(outcome == "territ_loss_concern"){
      basespec <-c("age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc",
                   "male", 
                   "national_language",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01",
                   "cultid",
                   "pplid",
                   "statid",
                   "town"
                   )
      if(country == "Hungary"){
        polit <- c("like_jobbik", 
                   "like_mszp",
                   "like_fidesz")
      } else if(country == "Romania"){
        polit <- c("like_prm",
                   "like_psd",
                   "like_pmp",
                   "like_udmr",
                   "like_usr",
                   "like_pnl")
      }else if(country == "Turkey"){
        polit <- c("like_akp", # Turkey
                   "like_mhp",
                   "like_hdp",
                   "like_chp",
                   "like_iyi")
      } else if(country == "Germany"){
        polit <- c("like_cducsu", 
                   "like_spd",
                   "like_afd")
      } else{
        polit <- c("populist_like_all_onep")
      }
      basespec <- c(basespec,polit)
  }

  f <- as.formula(
    paste(outcome, 
          ifelse(fe,paste(c(basespec,geo), collapse = " + "),paste(basespec, collapse = " + ")), 
          sep = " ~ "))
  return(f)
  
}
getmodels <- function(outcomes,countries,binary=F,controls=T,interact = F, obs=F,diaspora=F,fe=T,scale=F,quadratic=F){
  models <- list()
  if(binary){
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hung %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else if(countries[i] == "Romania"){
        models[[i]] <- roma %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turk %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else if(countries[i] == "Germany"){
        models[[i]] <- germ %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else{
        models[[i]] <- svydat %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      }
    }
  } else if (outcomes[[1]] == "territ_loss_concern"){
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hungall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Romania"){
        models[[i]] <- romaall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turkall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Germany"){
        models[[i]] <- germall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else{
        models[[i]] <- svydatall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }
    }
  } else{
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hung %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Romania"){
        models[[i]] <- roma %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turk %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Germany"){
        models[[i]] <- germ %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else{
        models[[i]] <- svydat %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }
    }
  }
  
  return(models)
}

#Created filtered survey weighted versions of data
svydat <-dat %>% filter(national_language==1) %>% svydesign(~1, strata = ~Country, weights=~wt,data=.)
svydatimpute <- subset(svydesign(~1, strata = ~Country, weights=~wt,data=svyimpute),national_language==1)
svydatimputehuntur <- subset(svydatimpute,hun==1|tur==1)
svydat4impute <-  subset(svydatimpute,Version!=4)
dat4impute <- subset_datlist(svyimpute,subset=c(datlist_create(svyimpute)[[1]]$Version!=4),toclass="imputationList")
svydat4imputehuntur <- subset(svydatimputehuntur,Version!=4)
svydat4 <-  subset(svydat,Version!=4)
svydat4votedpp <-  subset(svydat4,populist_voted==1)
svydat4novotepp <-  subset(svydat4,populist_voted==0)
svydatvotedpp <-  subset(svydat,populist_voted==1)
svydatnovotepp <-  subset(svydat,populist_voted==0)
svydatall <-dat %>% svydesign(~1, strata = ~Country, weights=~wt,data=.)

###################Summary statistics -- Table F.1##########################
t <- tbl_svysummary(svydat, include= c(
  "populist_like_all_max",
  "populist_like_in_power",
  "populist_like_not_in_power",
  "populist_voted",
  "hist_greatestextent",
  "territ_loss_concern",
  "age", 
  "educ_univ_complete", 
  "educ_sec_univprep", 
  "educ_sec_techvoc", 
  "male", 
  "city",
  "town",
  "covid_civlib_restrict_ind",
  "econ_intervention_ind",
  "pro_global_ind",
  "democ_ind",
  "natid_ind",
  "neigh_lgbt",
  "relig_import01",
  "pol_int01"),
  statistic=list(all_continuous() ~ "& {mean} & {sd} & {min} & {max} \\\\", all_categorical() ~ "0.{p}"), missing = "no",
  type = list(everything() ~ "continuous"))
#################################################################################
#Subset to indivdual countries
hung <- dat %>% filter(hun==1,!is.na(County_Hung),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
roma <- dat %>% filter(rom==1,!is.na(County_Rom),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
turk <- dat %>% filter(tur==1,!is.na(Province_Turk),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
germ <- dat %>% filter(ger==1,!is.na(Bundesland_Germ),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
hungimpute <- subset(svydatimpute,hun==1)
romaimpute <- subset(svydatimpute,rom==1)
turkimpute <- subset(svydatimpute,tur==1)
germimpute <- subset(svydatimpute,ger==1)

#######################Figure 1 -- data distributions############################
dat %>% filter(hun==1) %>% ggplot(.,aes(territ_loss_concern)) + 
  geom_histogram(aes(y = (..count..)/sum(..count..)),bins=11) + 
  theme_classic() + 
  xlab("Concern For Territorial Loss") +
  ylab("Proportion")
ggsave("../Figures/hundist.png",width=1920,height=1080,units="px")
dat %>% filter(rom==1) %>% ggplot(.,aes(territ_loss_concern)) + 
  geom_histogram(aes(y = (..count..)/sum(..count..)),bins=11) + 
  theme_classic() + 
  xlab("Concern For Territorial Loss") +
  ylab("Proportion")
ggsave("../Figures/romdist.png",width=1920,height=1080,units="px")
dat %>% filter(tur==1) %>% ggplot(.,aes(territ_loss_concern)) + 
  geom_histogram(aes(y = (..count..)/sum(..count..)),bins=11) + 
  theme_classic() + 
  xlab("Concern For Territorial Loss") +
  ylab("Proportion")
ggsave("../Figures/turdist.png",width=1920,height=1080,units="px")
dat %>% filter(ger==1) %>% ggplot(.,aes(territ_loss_concern)) + 
  geom_histogram(aes(y = (..count..)/sum(..count..)),bins=11) + 
  theme_classic() + 
  xlab("Concern For Territorial Loss") +
  ylab("Proportion")
ggsave("../Figures/gerdist.png",width=1920,height=1080,units="px")

#Versions without condition 4
hung4 <- dat4 %>% filter(hun==1,!is.na(County_Hung),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
roma4 <- dat4 %>% filter(rom==1,!is.na(County_Rom),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
turk4 <- dat4 %>% filter(tur==1,!is.na(Province_Turk),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
germ4 <- dat4 %>% filter(ger==1,!is.na(Bundesland_Germ),national_language==1) %>% svydesign(~1,weights=~wt,data = .)

hung4impute <- subset(svydat4impute,hun==1)
roma4impute <- subset(svydat4impute,rom==1)
turk4impute <- subset(svydat4impute,tur==1)
germ4impute <- subset(svydat4impute,ger==1)

hungall <- dat %>% filter(hun==1,!is.na(County_Hung)) %>% svydesign(~1,weights=~wt,data = .)
romaall <- dat %>% filter(rom==1,!is.na(County_Rom)) %>% svydesign(~1,weights=~wt,data = .)
turkall <- dat %>% filter(tur==1,!is.na(Province_Turk)) %>% svydesign(~1,weights=~wt,data = .)
germall <- dat %>% filter(ger==1,!is.na(Bundesland_Germ)) %>% svydesign(~1,weights=~wt,data = .)

######################################################################################
###############Party evals as function of concern -- Table D.1########################
######################################################################################
models <- list(
  with(hung4impute,svyglm(makeform("Hungary","like_jobbik",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(hung4impute,svyglm(makeform("Hungary","like_fidesz",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(hung4impute,svyglm(makeform("Hungary","like_mszp",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(roma4impute,svyglm(makeform("Romania","like_psd",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(roma4impute,svyglm(makeform("Romania","like_usr",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(roma4impute,svyglm(makeform("Romania","like_usr",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(roma4impute,svyglm(makeform("Romania","like_pnl",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(turk4impute,svyglm(makeform("Turkey","like_akp",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(turk4impute,svyglm(makeform("Turkey","like_chp",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(turk4impute,svyglm(makeform("Turkey","like_iyi",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(germ4impute,svyglm(makeform("Germany","like_afd",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(germ4impute,svyglm(makeform("Germany","like_spd",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T))),
  with(germ4impute,svyglm(makeform("Germany","like_cducsu",interact = T,
                                            controls=T,
                                            fe=T,
                                            scale = T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest",
                             "Prime : Concern"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Main Models########################
#########################################################################

# Table C.2
models <- list(
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",
                  obs=F,
                  controls=F,
                  fe=F))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",
                  controls=F,
                  fe=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",
                  controls=T,
                  fe=T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",
                  obs=F,
                  controls=F,
                  fe=F))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",
                  controls=F,
                  fe=T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",
                  controls=T,
                  fe=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",
                  obs=F,
                  controls=F,
                  fe=F))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",
                  controls=F,
                  fe=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",
                  controls=T,
                  fe=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")

texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Prime",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)
#Without survey weights -- Table C.3
models <- list(
  with(dat4impute,lm(populist_like_all_max ~ territorial)),
  with(dat4impute,lm(populist_like_all_max ~ territorial + Country)),
  with(dat4impute,lm(populist_like_all_max ~ territorial + age + educ_univ_complete + 
                                  educ_sec_univprep + educ_sec_techvoc + male + city + town + 
                                  covid_civlib_restrict_ind + econ_intervention_ind + pro_global_ind + 
                                  democ_ind + natid_ind + neigh_lgbt + relig_import01 + pol_int01 + 
                                  Country)),
  with(dat4impute,lm(populist_like_in_power ~ territorial)),
  with(dat4impute,lm(populist_like_in_power ~ territorial + Country)),
  with(dat4impute,lm(populist_like_in_power ~ territorial + age + educ_univ_complete + 
                                  educ_sec_univprep + educ_sec_techvoc + male + city + town + 
                                  covid_civlib_restrict_ind + econ_intervention_ind + pro_global_ind + 
                                  democ_ind + natid_ind + neigh_lgbt + relig_import01 + pol_int01 + 
                                  Country)),
  with(dat4impute,lm(populist_like_not_in_power ~ territorial)),
  with(dat4impute,lm(populist_like_not_in_power ~ territorial + Country)),
  with(dat4impute,lm(populist_like_not_in_power ~ territorial + age + educ_univ_complete + 
                                  educ_sec_univprep + educ_sec_techvoc + male + city + town + 
                                  covid_civlib_restrict_ind + econ_intervention_ind + pro_global_ind + 
                                  democ_ind + natid_ind + neigh_lgbt + relig_import01 + pol_int01 + 
                                  Country))
  )
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,function(y){summary(y)$r.squared}))})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Prime",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#Table C.1
models <- list(
  with(svydatimpute,svyglm(makeform("All","populist_like_all_max",
                  obs=T,
                  controls=F,
                  fe=F,
                  scale=T))),
  with(svydatimpute,svyglm(makeform("All","populist_like_all_max",
                  obs=T,
                  controls=F,
                  fe=T,
                  scale=T))),
  with(svydatimpute,svyglm(makeform("All","populist_like_all_max",
                  obs=T,
                  controls=T,
                  fe=T,
                  scale=T))),
  with(svydatimputehuntur,svyglm(makeform("All","populist_like_in_power",
                  obs=T,
                  controls=F,
                  fe=F,
                  scale=T))),
  with(svydatimputehuntur,svyglm(makeform("All","populist_like_in_power",
                  obs=T,
                  controls=F,
                  fe=T,
                  scale=T))),
  with(svydatimputehuntur,svyglm(makeform("All","populist_like_in_power",
                  obs=T,
                  controls=T,
                  fe=T,
                  scale=T))),
  with(svydatimpute,svyglm(makeform("All","populist_like_not_in_power",
                  obs=T,
                  controls=F,
                  fe=F,
                  scale=T))),
  with(svydatimpute,svyglm(makeform("All","populist_like_not_in_power",
                  obs=T,
                  controls=F,
                  fe=T,
                  scale=T))),
  with(svydatimpute,svyglm(makeform("All","populist_like_not_in_power",
                  obs=T,
                  controls=T,
                  fe=T,
                  scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")

texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
       )


#########################################################################
###############Interaction Models -- Table C.4 ##########################
#########################################################################


models <- list(
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=F,
                  controls=F,
                  fe=F,
                  scale=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=T,
                  controls=F,
                  fe=F,
                  scale=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=F,
                  controls=F,
                  scale=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=T,
                  controls=F,
                  scale=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=F,
                  controls=T,
                  scale=T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,
                  quadratic=T,
                  controls=T,
                  scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Prime : Concern",
                             "Territorial Loss Concern$^2$",
                             "Prime : Concern$^2$",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

####################################################################################
###############Interaction Models (in power) -- Table C.5 ##########################
####################################################################################


models <- list(
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=F,
                  controls=F,
                  fe=F,
                  scale = T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=T,
                  controls=F,
                  fe=F,
                  scale = T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=F,
                  controls=F,
                  scale = T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=T,
                  controls=F,
                  scale = T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=F,
                  controls=T,
                  scale = T))),
  with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,
                  quadratic=T,
                  controls=T,
                  scale = T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Prime : Concern",
                             "Territorial Loss Concern$^2$",
                             "Prime : Concern$^2$",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)
####################################################################################
###############Interaction Models (in opposition) -- Table C.6 #####################
####################################################################################


models <- list(
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=F,
                  controls=F,
                  fe=F,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=T,
                  controls=F,
                  fe=F,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=F,
                  controls=F,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=T,
                  controls=F,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=F,
                  controls=T,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,
                  quadratic=T,
                  controls=T,
                  scale = T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Prime : Concern",
                             "Territorial Loss Concern$^2$",
                             "Prime : Concern$^2$",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Historical Period Evaluations Pooled -- Table C.7#########
#########################################################################

models <- list(
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  controls=F,
                  fe=F),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  obs=T,
                  controls=F,
                  fe=F,
                  scale = T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = T,
                  controls=F,
                  fe=F,
                  scale = T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  controls=F,
                  fe=T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  obs=T,
                  controls=F,
                  fe=T,
                  scale = T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = T,
                  controls=F,
                  fe=T,
                  scale = T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  controls=T,
                  fe=T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = F,
                  obs=T,
                  controls=T,
                  fe=T,
                  scale = T),family=quasibinomial())),
  with(svydatimpute,svyglm(makeform("All","hist_greatestextent",interact = T,
                  controls=T,
                  fe=T,
                  scale = T),family=quasibinomial()))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("Pseudo R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Prime",
                             "Territorial Loss Concern",
                             "Prime : Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Historical Period Evaluations by country -- Table C.8#####
#########################################################################
models <- list(
  
  with(hungimpute,svyglm(makeform("Hungary","hunhist_pre18",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(hungimpute,svyglm(makeform("Hungary","hunhist_pre18",interact = T,
                                              controls=T,
                                              fe=T,
                                              scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","romhist_1938",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","romhist_1938",interact = T,
                                              controls=T,
                                              fe=T,
                                              scale = T),family=quasibinomial())),
  with(turkimpute,svyglm(makeform("Turkey","turhist_pre1922",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(turkimpute,svyglm(makeform("Turkey","turhist_pre1922",interact = T,
                                              controls=T,
                                              fe=T,
                                              scale = T),family=quasibinomial())),
  with(germimpute,svyglm(makeform("Germany","gerhist_reich",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(germimpute,svyglm(makeform("Germany","gerhist_reich",interact = T,
                                              controls=T,
                                              fe=T,
                                              scale = T),family=quasibinomial()))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("Pseudo R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest",
                             "Prime : Concern"),
       stars=c(0.01,0.05,0.1)
)


#########################################################################
########################## Vote Choice -- Table D.2######################
#########################################################################
models <- list(
  with(hungimpute,svyglm(makeform("Hungary","voted_jobbik",
                                            obs=T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(hungimpute,svyglm(makeform("Hungary","voted_fidesz",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(hungimpute,svyglm(makeform("Hungary","voted_nonntl_hun",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","voted_dancila",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","voted_paleologu",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","voted_barna",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(romaimpute,svyglm(makeform("Romania","voted_iohannis",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(turkimpute,svyglm(makeform("Turkey","voted_akp",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(turkimpute,svyglm(makeform("Turkey","voted_chp",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(turkimpute,svyglm(makeform("Turkey","voted_iyi",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(germimpute,svyglm(makeform("Germany","voted_afd",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(germimpute,svyglm(makeform("Germany","voted_spd",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial())),
  with(germimpute,svyglm(makeform("Germany","voted_cdu",obs = T,
                                            controls=T,
                                            fe=T,
                                            scale = T),family=quasibinomial()))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

############################################################################################
########################## Coefficient Plots for pooled results ############################
############################################################################################

#Helper function for marginal plots
miplot <- function(models,model.names,coefs,ci_level,colors,...){
  sink("file")
  estimates <- unlist(lapply(models,function(x){return(summary(x)$results[match(coefs,names(x$coefficients))])}))
  ses <- unlist(lapply(models,function(x){return(summary(x)$se[match(coefs,names(x$coefficients))])}))
  sink()
  upper <- estimates +qnorm(1-(1-ci_level[1])/2)*ses
  lower <- estimates -qnorm(1-(1-ci_level[1])/2)*ses
  mnames <- unlist(sapply(model.names,function(x){rep(x,length(coefs))},simplify=F))
  if(length(names(coefs))>0){outnames <- names(coefs)} else{outnames <- coefs}
  
  data <- tibble(estimate=estimates,std.error=ses,conf.high=upper,conf.low=lower,term = rep(outnames,length(models)),model=mnames)
  if(length(ci_level>1)){
    z <- qnorm(1-(1-ci_level[2])/2)
    data <- data %>% mutate(clow2 = (estimate - z*std.error),chigh2 = (estimate + z*std.error)) 
  }
  p<- dwplot(data,...)+ labs(color='Model') + 
    theme_bw() +
    scale_color_manual(values=colors) + 
    geom_vline(xintercept =0)
  if(length(ci_level>1)){
    jitter <- ggplot_build(p)$data[[1]]$y
    p <- p + geom_segment(aes(x=clow2,y=jitter,xend=chigh2,
                              yend=jitter,col=model),size=0.5)
  }
  return(p)
}

#Figure 4(a)
 m1 <- with(svydat4impute,svyglm(makeform("All","populist_like_all_max",controls=F,fe=F)))

 m2 <- with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",controls=F,fe=F)))

 m3 <- with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",controls=F,fe=F)))

p <- miplot(lapply(list(m1,m2,m3),MIcombine),coefs=c("Loss Prime" = "territorial1"),
                model.names =c("All","In Power", "In Opposition"),
                ci_level=c(0.9,0.95),
                colors = c("#fde725","#21918c","#440154"),whisker_args=list(size=1))
p
ggsave("../Figures/Experimental Effects.png",width=1920,height=1080,units="px")

#Figure 4(b)
m1 <- with(svydatimpute,svyglm(makeform("All","populist_like_all_max",obs=T)))
m2 <- with(svydatimputehuntur,svyglm(makeform("All","populist_like_in_power",obs=T)))
m3 <- with(svydatimpute,svyglm(makeform("All","populist_like_not_in_power",obs=T)))
p <- miplot(lapply(list(m1,m2,m3),MIcombine),coefs=c("Loss Concern" = "territ_loss_concern"),
            model.names =c("All","In Power", "In Opposition"),
            ci_level=c(0.9,0.95),
            colors = c("#fde725","#21918c","#440154"),whisker_args=list(size=1))
p
ggsave("../Figures/Observational Effects.png",width=1920,height=1080,units="px")



m1 <- MIcombine(with(svydat4impute,svyglm(makeform("All","populist_like_all_max",interact = T,scale=F,quadratic=T,controls=T))))

 m2 <- MIcombine(with(svydat4imputehuntur,svyglm(makeform("All","populist_like_in_power",interact = T,scale=F,quadratic=T,controls=T))))

 m3 <- MIcombine(with(svydat4impute,svyglm(makeform("All","populist_like_not_in_power",interact = T,scale=F,quadratic=T,controls=T))))

#Figure 5
estimates1 <- meplot(m1,"territorial1","territ_loss_concern")

ggplot(estimates1,aes(x=x,y=est)) +
  geom_point() +
  geom_errorbar(aes(ymax=ci.upper,ymin=ci.lower)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  xlab("Concern for Territorial Loss") +
  ylab("Marginal Effect of Prime") +
  ylim(c(-35,15))
ggsave("../Figures/marginaleffectsall.png",width=1920,height=1080,units="px")
#Figure 6(a)
estimates2 <- meplot(m2,"territorial1","territ_loss_concern")

ggplot(estimates2,aes(x=x,y=est)) +
  geom_point() +
  geom_errorbar(aes(ymax=ci.upper,ymin=ci.lower)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  xlab("Concern for Territorial Loss") +
  ylab("Marginal Effect of Prime")
ggsave("../Figures/marginaleffectsinpower.png",width=1920,height=1080,units="px")
#Figure 6(b)
estimates3 <- meplot(m3,"territorial1","territ_loss_concern")

ggplot(estimates3,aes(x=x,y=est)) +
  geom_point() +
  geom_errorbar(aes(ymax=ci.upper,ymin=ci.lower)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  xlab("Concern for Territorial Loss") +
  ylab("Marginal Effect of Prime")
ggsave("../Figures/marginaleffectsnotinpower.png",width=1920,height=1080,units="px")

#Figure E.1
estimates <- bind_rows(list(estimates1,estimates2,estimates3),.id="model") %>% 
  mutate(Outcome = case_match(model,
    "1" ~ "All",
    "2" ~ "In Power",
    "3" ~ "In Opposition"
  ))

ggplot(estimates,aes(x=x,y=est,col=Outcome,shape=Outcome)) +
  geom_point(alpha=0.5) +
  geom_errorbar(alpha=0.3,aes(ymax=ci.upper,ymin=ci.lower)) +
  theme_bw() +
  geom_hline(yintercept = 0) +
  xlab("Concern for Territorial Loss") +
  ylab("Marginal Effect of Prime") +
  scale_color_viridis(discrete = T)
ggsave("../Figures/marginaleffectscombined.png",width=1920,height=1080,units="px")

#Figure 7
m <- svyglm(makeform("All","hist_greatestextent",interact = T,scale=F,quadratic = F),design=svydat,family = quasibinomial())

p <- plot_model(m,type = "pred",terms=c("Loss Concern" = "territ_loss_concern","Loss Prime" = "territorial"))
p+ labs(color='Loss Prime Status') + 
  theme_bw() +
  scale_color_manual(values=c("#440154","#21918c")) +
  scale_fill_manual(values=c("#440154","#21918c")) +
  xlab("Concern for Territorial Loss") +
  ylab("Predicted Probability") + 
  ggtitle("")
ggsave("../Figures/Historical Period Predicted.png",width=1920,height=1080,units="px")

################################################################################
########################## AUR Panel Analyses -- Tables 4 and C.9###############
################################################################################

p13[p13==99] <- NA
p13[p13==98] <- NA

p13$wt <- p13$pred/1000
m1 <- p13  %>% mutate(territ_concern = territ_concern/100,w1coveval=w1coveval/4) %>% glm(vote_aur ~ territ_concern + 
                                                                                           bessarabia + 
                                                                                           covid_civlib_restrict_ind + 
                                                                                           age + 
                                                                                           female + 
                                                                                           univ_edu + 
                                                                                           city + 
                                                                                           econ_intervention_ind +
                                                                                           natid_ind + 
                                                                                           neigh_lgbt + 
                                                                                           relig_import01 +
                                                                                           pol_int01 +
                                                                                           as.factor(County_Rom_w1),data=.,family=quasibinomial(),weights=wt)
m2 <- p13  %>% mutate(territ_concern = territ_concern/100,w1coveval=w1coveval/4) %>% lm(pty_close_AUR ~ territ_concern + 
                                                                                          bessarabia + 
                                                                                          covid_civlib_restrict_ind + 
                                                                                          age + 
                                                                                          female + 
                                                                                          univ_edu + 
                                                                                          city + 
                                                                                          econ_intervention_ind +
                                                                                          natid_ind + 
                                                                                          neigh_lgbt + 
                                                                                          relig_import01 +
                                                                                          pol_int01 +
                                                                                          as.factor(County_Rom_w1),data=.,weights=wt)

m3 <- p34 %>% mutate(territ_loss_concern_1 = territ_loss_concern_1/100,territ_loss_concern_4 = territ_loss_concern_4/100) %>% lm(territ_loss_concern_4 ~ vote_aur+  territ_loss_concern_1 + age_w3 + female_w3 + univ_edu_w3 + city_w3,data=.)

stargazer(list(m1,m2,m3),
          style="apsr",
          omit = "County",
          notes = c("All models include regional fixed effects and inverse probability attrition weights."), 
          notes.append = T)

#Table C.8
round(svymean(~factor(hist),hung) *100,1)
round(svymean(~factor(hist),roma) *100,1)
round(svymean(~factor(hist),turk) *100,1)
round(svymean(~factor(hist),germ) *100,1)
round(svymean(~sudetenland,germ) *100,1)

